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Concepts underlying the Enskog kinetic theory of hard-spheres are applied to include short-range 
correlation effects in a model for transport coefficients of strongly coupled plasmas. The approach 
is based on an extension of the effective potential transport theory [S. D. Baalrud and J. Daligault, 
Phys. Rev. Lett. 110 , 235001 (2013)] to include an exclusion radius surrounding individual charged 
particles that is associated with Coulomb repulsion. This is obtained by analogy with the finite size 
of hard spheres in Enskog’s theory. Predictions for the self-diffusion and shear viscosity coefficients of 
the one-component plasma are tested against molecular dynamics simulations. The theory is found 
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I. INTRODUCTION 

Plasmas found in several modern research areas can 
span a broad range of Coulomb coupling strengths. Ex¬ 
amples include dusty, non-neutral, and ultracold plasma 
experiments, as well as dense plasmas found in inertial 
confinement fusion, high-intensity laser-matter interac¬ 
tion experiments and dense astrophysical objects. A re¬ 
search need common to all of these areas is a transport 
theory that is versatile enough to cover a broad range 
of coupling strengths, and can be evaluated efficiently 
enough to be incorporated into the fluid codes used to 
simulate these systems. To date, no systematic theory 
is available to do this, yet a workable theory may be ob¬ 
tained through ad hoc extensions of theories of simplified 
systems. 

We recently proposed an approach based on an effec¬ 
tive interaction potential in an effort to extend conven¬ 
tional plasma transport theory into the strongly coupled 
regime [1]. This is a physically-motivated approach based 
on a Boltzmann-likc binary collision operator, but where 
many-body correlation effects are modeled through an ef¬ 
fective interaction potential. By comparing with molec¬ 
ular dynamics (MD) simulations of diffusion and shear 
viscosity of the one-component plasma (OCP), this was 
shown to be successful at extending the binary collision 
approach well into the strongly coupled regime. How¬ 
ever, one persistent feature was a 30-40% underestima¬ 
tion of the collision rate in the range 1 < T < 30, 
where T = e 2 /(aksT) is the Coulomb coupling parame¬ 
ter, a = ( 3747 m ) 1 / 3 is the Wigner-Seitz radius and T the 
temperature Hi- The theory begins to break down at 
larger coupling strengths (T > 30) where direct interac¬ 
tion contributions that are not included in a Boltzmann- 
like treatment begin to dominate. In this paper, we at¬ 
tempt to improve the theory by invoking ideas underly¬ 
ing the kinetic theory proposed by Enskog to extend the 
kinetic theory of Boltzmann to dense gases [5J [6j. The 
Enskog kinetic equation was developed for hard-spheres 


only and involves the introduction of corrections to the 
Boltzmann equation that account for the finite particle 
size. Although no one has yet succeeded in deriving a 
similar equation for continous potentials, the concepts 
underlying Enskog’s theory provide a valuable aid in un¬ 
derstanding the physical origin of the underestimation of 
the collision rate, as well as a source of ideas to improve 
the effective potential theory. 

Two basic assumptions limit the Boltzmann equa¬ 
tion to dilute systems: (i) only pairs of particles col¬ 
lide simultaneously (i.e., binary collisions), and (ii) 
the “Stosszahlansatz” or “molecular-chaos” assumption. 
The effective potential theory relaxes assumption (i) 
somewhat by treating binary scatterers as interacting 
through the potential of mean force, rather than the bare 
Coulomb potential. The potential of mean force is ob¬ 
tained by taking the two scattering particles at fixed po¬ 
sitions and averaging over the positions of all other par¬ 
ticles \T \. It is related to the pair distribution function by 
7>(r)/fcsT = — ln[g(r)], and includes many-body effects 
of the background including screening and correlations. 
The present work shows that further extension may be 
realized by addressing assumption (ii) through Enskog’s 
equation. 

Boltzmann’s molecular chaos approximation assumes 
that particles are uncorrelated prior to a collision, and 
even ignores the difference in the positions of the two 
colliding particles by setting n = r 2 = r. Enskog’s the¬ 
ory relaxes these assumptions by modeling molecules as 
hard spheres of finite diameter cr; see Fig. |TJ This treats 
the fact that particle centers cannot be closer than their 
physical diameter (an exclusion radius), which introduces 
both an aspect of correlation in the distribution of initial 
positions of scattering particles, as well as nonlocal as¬ 
pects related to the scale at which momentum and energy 
transfer occurs in a collision. These give rise to qualita¬ 
tively new features in the resulting fluid theory, including 
a non-ideal equation of state, and potential energy con¬ 
tributions to transport coefficients such as shear viscosity 
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and thermal conductivity [Bi. 

Enskog’s equation provides great insight, but real gases 
are not comprised of hard spheres. The hard sphere ap¬ 
proximation is central in a systematic derivation of En¬ 
skog’s equation [2j, but successful descriptions of realistic 
fluids have been realized through ad hoc identification of 
appropriate effective particle diameters [3j. The present 
work seeks a parallel approach for plasmas. At the out¬ 
set, it is not immediately apparent that this is possible 
since the Coulomb potential is “soft” in the sense that it 
falls off gradually (oc 1/r) from particle centers. A crit¬ 
ical aspect of the present work is the effective potential 
concept: colliding pairs interact through the potential of 
mean force rather than the bare Coulomb potential. At 
strong coupling many-body correlations cause the poten¬ 
tial of mean force to sharply fall-off at a fixed distance 
from the particle center, forming a “Coulomb hole”; see 
Fig- [?} We show that this can be associated with an effec¬ 
tive particle radius, opening a path for applying Enskog’s 
theory. 

Two parameters arise in Enskog’s modification to the 
dilute gas transport coefficients: An increased collision 
frequency (x) and the hard sphere packing fraction ( fj ). 
There is no systematic way to map these hard sphere 
parameters onto plasmas. We first determine them by 
equating MD simulations of the OCP self diffusion and 
viscosity coefficients with Enskog’s expressions, taking 
the effective potential theory as the dilute gas compo¬ 
nent. Essentially no modification from the dilute gas 
effective potential theory is found for T < 1. The 
main modification in the region 1 ^ T < 30 is an in¬ 
creased collision frequency of approximately 30 — 40% 
(x — 1.3 — 1.4). Several attempts to model these param¬ 
eters directly from the pair distribution function g(r ) are 
described. Although attempts at a systematic derivation 
were largely unable to quantitatively predict the Enskog 
parameters, a heuristic method based on associating the 
particle size with the Coulomb hole in g(r ) gives strong 
support for the notion of an increased collision frequency 
associated with an exclusion radius. 

The practical advantage of this approach is that it pro- 



FIG. 1. Geometry of two colliding hard spheres of diameter 
a at the point of contact. 



FIG. 2. (color online) (a) Pair distribution function for hard 
spheres calculated from the Percus-Yevick approximation, (b) 
Pair distribution function the OCP calculated from the HNC 
approximation (lines) and MD simulations (circles). 


vides a computationally efficient and accurate way to de¬ 
termine transport coefficients in complicated systems. In 
comparison, MD simulations are known to provide highly 
accurate results, but are sufficiently computationally ex¬ 
pensive that the prospect of building tables of transport 
coefficients over a large range of densities, temperatures, 
and mixture concentrations quickly becomes impracti¬ 
cal. In this paper, the theory is validated against MD 
simulations of simple OCP and Yukawa OCP systems. 
Applications to more realistic systems such as mixtures 
and dense plasmas will be published subsequently. 

The remainder of this paper is organized as follows. 
Section [n] briefly reviews the molecular chaos assump¬ 
tion in Boltzmann’s theory, and Enskog’s kinetic theory 
for the dense hard sphere gas. The simplification pro¬ 
vided by the hard-sphere approximation is that binary 
collisions take place instantaneously at known relative 
positions. This is not the case with continuos interparti¬ 
cle positions and the difficulties that arise in a systematic 
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FIG. 3. Geometry of a binary collision between point parti¬ 
cles in the center of mass frame. 

generalization of Enskog’s kinetic theory to treat par¬ 
ticles interacting through continuous potentials are dis¬ 
cussed in Sec. IIII1 This section also introduces a modified 
Enskog transport theory based on an effective particle 
diameter. Section |IV| provides a calculation of the En¬ 
skog parameters from MD simulations, and compares the 
results with attempts to determines these directly from 
g{r). Section |V| provides a comparison between the pre¬ 
vious effective potential theory calculations, calculations 
based on the modified Enskog theory, and MD simulation 
data for self-diffusion and shear viscosity of the OCP. Im¬ 
plications of these results are discussed in Sec. m 

II. ENSKOG THEORY 

A. Boltzmann’s Stosszahlansatz 

In the traditional derivation of the Boltzmann equa¬ 
tion, one calculates the expected number of binary col¬ 
lisions dN experienced during a small time interval dt 
by particles ‘1’ located in the phase-space volume el¬ 
ement d 3 rid 3 vi centered about the phase-space point 
(ri,vi) 0. For such a collision to occur, the collid¬ 
ing partners ‘2’ with velocity in the range dv 2 around 
v 2 must initially lie within a cylinder, the so-called col¬ 
lision cylinder, having an area bdbd(j) and generator 
— (vi — v 2 )df. Its volume is dV = |vi — \r 2 \b db dc/) dt; 
see Fig. [3] Here, b is the impact parameter and <f> is 
the azimuthal angle. In general, the expected number 
of binary collisions dN is the integral over the collision 
cylinder of the total number of pairs of collision partners 

dN= f 2 (r 1 ,v 1 ;r 2 ,v 2 ;t)dr 2 dv 2 dr 1 dv 1 , (1) 

Jr 2 &dV 

where / 2 is the phase-space pair distribution function. 
The Boltzmann equation additionally relies on the molec¬ 
ular chaos assumption (a.k.a. Stosszahlansatz) according 
to which (i) colliding particles are entirely uncorrelated 


in position before colliding, i.e. / 2 (iq, Vi; r 2 , v 2 ; t) ft 
/i(ri,vi;t)/i(r 2 , v 2 ;t), and (ii) the difference in position 
between two molecules is ignored and the distribution 
functions /1 are evaluated at the same point iq = r 2 = r 
in space. The molecular chaos approximation can then 
be written 

/ 2 (ri,vi;r 2 ,v 2 ;f) ft /i(r,Vi,t)/i(r, v 2 ,t) , (2) 

where f± is the single-particle distribution function. 
Combining this with Eq. 0 yields 

dN ft /i(r, vi,t)/i(r, v 2 ,t)|vi — v 2 \b db d(f> dtd 3 v 2 d 3 rd' i v 1 

The molecular chaos assumption restricts the validity of 
the Boltzmann equation to low densities, or in the con¬ 
text of the present paper, to small coupling strengths. 

In general, particle positions are correlated and corre¬ 
lations increase with density or coupling strength. For 
instance, in thermal equilibrium 

/ 2 (ri,vi;r 2 ,v 2 ;i) = g{ iq - r 2 )/ MB (vi)/ M B(v 2 ) (3) 

where /mb is the Maxwell-Boltzmann distribution and g 
is the pair distribution function. Physically, ng(r) rep¬ 
resents the radial density around a given particle. In 
essence, molecular chaos assumes that g(r) = 1 for all 
distances r. Figure [2j;> shows the pair distribution func¬ 
tion of the OCP at several coupling strengths over a range 
of distances on the order of the interparticle spacing. For 
weak coupling (r < 0.1), g{r) is indeed very close to one 
over the range relevant to colliding particles. However, 
as the coupling strength increases, g(r ) transitions from 
0 to 1 at a radial distance on the order of the interparticle 
spacing. It stays nearly equal to zero over an increasingly 
large range as T increases and then oscillates around one. 
The “hole” in g(r) at small r can be interpreted as an 
exclusion volume around each particle, which originates 
from the strong repulsion at small distances that the vast 
majority of particles can not penetrate due to their low 
kinetic energy. Figure [2[i illustrates that a similar exclu¬ 
sion zone is found for hard spheres, which naturally arises 
from the finite size of particles. In this case, the exclu¬ 
sion zone is simply related to the size of a particle and 
is independent of temperature and density; the value of 
g(r) depends on the physical conditions only beyond the 
particle diameter (r > a). On the contrary, with parti¬ 
cles interacting via a continuous potential, the exclusion 
zone depends on the physical conditions and cannot be 
delineated as clearly. 

B. Enskog’s Equation 

In 1922 Enskog was able to successfully abandon the 
molecular chaos approximation for hard sphere gases [3] . 
He replaced it by an ansatz that includes effects of spa¬ 
tial correlations by invoking the following arguments: (i) 
Because of their finite size, the centers of colliding hard 
spheres are not at the same position at contact. Rather, 
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they are separated by a distance equal to the particle di¬ 
ameter. Thus, /i(r, vi,t)/i(r,v 2 ,f) should be replaced 
by 

/i(r, vi,i)/i(r — crk, v 2 ,t) (4) 

in the molecular chaos approximation, Eq. ©• Here, k is 
the unit vector joining the hard-sphere centers; see Fig. [l] 
(ii) In addition, the finite particle size reduces the avail¬ 
able volume that particles can occupy, which increases 
the probability of a collision. Enskog accounted for this 
by multiplying Eq. © by a factor y(r — ^crk) evaluated 
at the point of contact, i.e. 

X( r ~ ^CTk)/i(r,Vi,f)/i(r-f7k,v 2 ,t). (5) 

Enskog derived his kinetic equation from similar 
heuristic arguments as used by Boltzmann, but replac¬ 
ing the molecular chaos approximation by Eq. © and 
considering exclusively hard spheres. The result was 
df /dt = Ce(/, /) with the collision operator 

Ce(/, /) = Jd 3 v'd 2 ka 2 u ■ k (6) 

x( r +^k)/(r)/( r +o-k) - x(r- *erk)/(r)/(r - ak) 

Here, u = vi — v 2 , and / denotes distribution func¬ 
tions evaluated at post-collision velocities v = v + Av, 
whereas / denotes distribution functions evaluated at 
pre-collision velocities v. In order for the resulting trans¬ 
port equations to be consistent with the equation of state, 
the factor y is taken equal to the equilibrium pair dis¬ 
tribution function % = g(o/2) at the point of contact 
|(ri + r 2 ) = f k between the two colliding particles. 

Equation (|6| differs from the Boltzmann collision op¬ 
erator in the following ways: (a) The factor y is absent 
in the Boltzmann collision operator, (b) The distribution 
functions are evaluated at non-local spatial locations in 
Eq. © , whereas all are evaluated at the local position r in 
the Boltzmann equation, (c) The kernel associated with 
the scattering probability in Eq. ©, d 2 ka 2 u k is partic¬ 
ular to hard spheres, whereas in the Boltzmann equation 
this is replaced by dQcr' u where a' is a differential scatter¬ 
ing cross section associated with an unspecified potential 
v(r). 

Enskog then derived fluid transport equations by ap¬ 
plying the Chapman-Enskog expansion to Eq. © . In 
comparison to the transport coefficients derived from the 
Boltzmann equation, aspect (a) leads to an increased col¬ 
lision frequency in the Enskog theory. This is manifest 
in the self-diffusion coefficient 

D = D 0 / x (7) 

where D 0 is the dilute gas self-diffusion coefficient 
obtained from the Boltzmann equation. For hard 
spheres, the lowest order expression is D 0 ^ s = 
^(nmkBT) 1 / 2 /(8nmncr 2 ) [Bj. 


In addition to the increased collision frequency, the 
non-local aspect (b) gives rise to a non-ideal equation of 
state, and contribution from particle interactions in the 
viscosity and thermal conductivity transport coefficients. 
Here, we will be interested in the viscosity coefficient 

V=~[l + 0.8 bp X + 0.7614(6py) 2 ] (8) 

X 

where rj 0 is the dilute gas shear viscosity. For 
hard spheres, the lowest order term is r/ 0t \ ls = 
5(7r?nA:BT) 1 / 2 /(167rtj 2 ). Here, 

2 

bp = -Trna 3 (9) 

O 

is the co-volume of the molecules, which can also be ex¬ 
pressed in terms of the packing fraction fj = bp/4. The 
first term on the right of Eq. ([8]) is the kinetic term, 
which is the dilute gas result modified only by the in¬ 
creased collision probability factor y. The second two 
terms on the right arise due to the non-local aspects of 
the collision operator. These potential contributions are 
absent in the dilute gas theory. 

Evaluation of these coefficients requires some external 
determination of y. This is usually obtained by equating 
the equation of state implied by Eq. ©. p/(nk B T) = 
1 + bpx, with that obtained by another means. One op¬ 
tion is the Carnahan-Starling approximation [5J, which 
leads to yes = (1 — fj/ 2) / (1 — fj) 3 . Similarly, the Percus- 
Yevick approximation m gives y PY = {l+g/2)/(l-rj 2 ). 
Another approach is to equate Enskog’s equation of state 
with the virial coefficients of the thermodynamic equa¬ 
tion of state for hard spheres p/(nk B T ) = 1 + bp + 
0.6250 (bp) 2 + 0.2869 (bp) 3 + 0.115(6p) 4 + .... This im¬ 
plies [6j 

y = 1 + 0.6250fep + 0.2869(6p) 2 + 0.115(6p) 3 + ... (10) 

The accuracy of these various approximations have been 
compared in detail (e.g., see mm)- 

It is not obvious how to extend Enskog’s ansatz, 
Eq. ©, to continuous interparticle potentials. For in¬ 
stance, how does one define an effective particle size? 
The size of a hard sphere is an intrinsic property, but for 
continuous potentials the exclusion volume is a statisti¬ 
cal concept and the distance of closest approach between 
particles in a binary collision is momentum dependent. 
Nevertheless, it is reasonable to expect that the corre¬ 
lation hole will affect transport properties in a similar 
fashion as the finite particle size of hard-spheres. In fact, 
in his 1922 paper [5], Enskog himself gave indications as 
to how his theory might be applied to real systems. In 
the next section, we discuss challenges that arise when 
attempting a systematic derivation of a generalized En¬ 
skog equation for soft potentials, then we outline a phe¬ 
nomenological approach that enables progress on approx¬ 
imating fluid transport coefficients. 
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FIG. 4. (a) Depiction of momentum vectors in a binary col¬ 

lision of hard spheres, and (b) in a binary collision of particles 
interacting through a continuous potential. 


III. MODIFIED ENSKOG THEORY 
A. Difficulties of a Systematic Theory 

To our knowledge, no one has yet succeeded in deriv¬ 
ing an Enskog-like kinetic equation for continuous poten¬ 
tials either phenomenologically or from first-principles. 
In fact, it is only in 1961 that Sengers and Cohen |8j could 
justify the Enskog kinetic equation from the BBGKY hi¬ 
erarchy. They did so using two independent methods: 
the coarse-graining method of Bogolyubov and the time¬ 
smoothing approach of Kirkwood. This section provides 
some identification of the difficulties that arise when at¬ 
tempting to extend such methods to treat continuous po¬ 
tentials. To this end, we highlight important steps of 
Sengers and Cohen’s derivation, giving a simplified ver¬ 
sion to make our point more transparent and accessible. 
We highlight when and why extending the derivation to 
continuous potentials faces serious difficulties. 

Following the coarse-graining method that Kirkwood 
used to derive the Boltzmann equation in his 1946 pa¬ 
per m, the distribution functions of the BBGKY hier¬ 


archy are time-smoothed over an interval from 0 to r, 
where the time r is long compared to the time of a colli¬ 
sion and short compared to the time between collisions. 
For instance, for the single-particle distribution function, 

fi(x,t) = ~ [ fi{x,t + s)ds (11) 

T Jo 

where x = (r, p) is a point in phase-space. The latter 
satisfies the first BBGKY equation, 


^ +V ^ = r/ dS j dx 2 Ll 2 -f 2 ( Xl ’ X2 ’ t + s ) ~ C 

( 12 ) 

where L\ 2 = — 5m) interaction 

vertex and f 2 is the two-particle distribution function. 
The latter satisfies the second BBGKY equation. Thus, 
within the binary collision approximation underlying the 
effective potential theory, we can write 

h(xi, x 2 \ t + s) = e~ sHl2 f 2 {xi,x 2 ]t) (13) 


z z 

where %i 2 . = ^ + 4>(r 12 ),.} is the Liouville oper¬ 

ator for a two-particle system. In Eq. (131, e~ sUl2 prop¬ 
agates in time the dynamics of two particles over a time 
interval r. Within this binary collision approximation, 
the right-hand side of Eq. (121 becomes 



dx 2 L 12 e sUl2 f 2 (x 1 ,x 2 ;t). (14) 


For hard-spheres, all binary collisions happen instan¬ 
taneously and one can take the time average over an in¬ 
finitesimally small time interval r —> 0 + during which 
only one binary collision can occur, and one has 

C = lim - I ds [ dx 2 L 12 e~ sHl2 f 2 (xi,x 2 ',t). (15) 
T^0+ T J Q J 

As illustrated in Fig. [1J the particle position does not 
change in an instantaneous collision, so 


lim e sWl2 ri = iq, lim e sWl2 r 2 = r 2 (16) 

T-90+ t—>0+ 


whereas the momenta change discontinuously 


lim e ® Wl2 p 1 =p 1 , lim e sWl2 p9 = p 2 . 

T-X)+ T—>0+ 


(17) 


Noting that r 2 


ri + uk at contact, Eq. (15) reduces to 


C = lim - [ ds [ dx 2 L 12 f 2 {r 1 ,e sWl2 pi;r 2 ,e sHl2 p 2 ;i) 

T-m+ r Jq J 

= Jd p 2 / dku 2 |(u • k)| [/ 2 (n, pi; n + ak, p 2 ; t) - / 2 (n, p, n + uk, p 2 ; t) ] 


(18) 

(19) 


where the last equation is derived in |5]. Now, if in addition one assumes that f 2 can be approximated by 
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h{x\,x 2 ]t) « g{Vi - r 2 )fi{xi,t)f 1 (x 2 ,t) (as in the equi¬ 
librium relation), the latter equation reduces to the En- 
skog kinetic equation. 

The extension of the previous steps to continuous po¬ 
tentials leads to serious difficulties. As illustrated in 
Fig. B (i) Collisions are not instantaneous, therefore the 
smoothing time scale s is finite and one needs to perform 
the time averaging integral in C. (ii) The duration of a 
collision depends on the collision parameters, therefore 
the time s depends on collision parameters, (iii) Both 
positions and momenta change during a binary collision, 
which must be incorporated in Eqs. (16])-(T7|). Overcom¬ 
ing these difficulties is a substantial, and longstanding, 
challenge in kinetic theory. However, practical progress 
can be made through analogies between the hard-sphere 
system and real systems at the level of macroscopic trans¬ 
port equations. The remainder of this paper explores this 
track. 


B. Modified Enskog Theory For Plasmas 

Already in his seminal paper of 1922, Enskog suggested 
how to adjust his theory to real neutral gases in an ad-hoc 
manner [5j; this approach is usually called the modified 
Enskog theory. The dilute gas coefficients are computed 
from the usual Chapman-Enskog solution of the Boltz¬ 
mann equation using the interaction potential associated 
with the realistic system. In addition, he suggested that 
the value of bpx be calculated to reproduce the varia¬ 
tion of pressure with temperature (known as the thermal 
pressure) instead of simply the pressure as in the original 
Enskog theory. This was found to be accurate especially 
for simple gases at packing fractions corresponding to 
bp < 0.6 (for example, see [UES]) 

Our goal is analogous to that of the modified Enskog 
theories in that we seek appropriate analogies to relate 
key features of the realistic system to the hard sphere 
model. However, a couple of difficulties are immediately 
apparent when the realistic system is a plasma rather 
than a neutral fluid. One is that the density is no longer 
the appropriate expansion parameter to link the weakly 
coupled to the strongly coupled regimes, and it should 
be replaced by the Coulomb coupling parameter. More¬ 
over, traditional plasma theory based on a dilute gas¬ 
like Boltzmann equation diverges in the strongly cou¬ 
pled regime. We propose that the appropriate dilute gas 
theory is provided by the effective potential theory from 
mm- Like the neutral fluid dilute gas theory, this is 
also based on the Boltzmann kinetic equation. The dif¬ 
ference is that the interaction potential is the potential 
of mean force, rather than the bare particle potential. 
Making this association, the corresponding self-diffusion 
coefficient is 


D* 


D i 

aksT 


= —D 
X 


EP- 


( 20 ) 


where Z? EP is the diffusion coefficient from the ef¬ 


fective potential theory. At lowest order, D\ EP = 

^73/(r 5 / 2 5(M)). Here, E^ 1,1 ) is a generalized Coulomb 
logarithm calculated from the effective potential the¬ 
ory mm- Similarly, the shear viscosity coefficient is 


V = 


Vi 


%p 

X 


[l+0Mpx + 0.76U(bp X ) 2 }. (21) 


At lowest order, the shear viscosity coefficient ob¬ 
tained from the effective potential theory is EP = 
5 v / 7r/(3v / 3T 5 / 2 S( 2 ’ 2 )), where S* 2 - 2 ) is a generalized 
Coulomb logarithm from the effective potential theory 
as discussed in M- 

Another apparent difficulty is that the effective size of 
particles is associated with a statistical “Coulomb hole” 
generated by electrostatic repulsion, rather than a hard 
core. Thus, a model is required to quantify the effective 
diameter of particles, as well as to obtain an approxi¬ 
mation for x- Typical relationships between model and 
hard sphere equations of state are not viable for plasmas. 
For instance, the OCP pressure becomes highly negative 
at strong coupling due to the presence of the homoge¬ 
neous neutralizing background m- Since both bp and 
X are positive definite quantities, the equation of state 
does not provide a meaningful relationship between the 
two systems. In the next section, we first determine the 
parameters bp and x i n Eqs. (201 and (21) from MD sim¬ 


ulations of D* and rj*. The results are then compared 
with attempts to approximate these parameters directly 
from g(r). 


IV. EFFECTIVE SIZE OF COULOMB POINT 
PARTICLES 

A. Constraints From MD Simulations 


If the modified Enskog theory described by Eqs. (20) 


and (21) are assumed to provide an accurate model, MD 
simulations can be used to determine the two indepen¬ 
dent parameters bp and ;\;. The enhanced collision prob¬ 
ability factor from Eq. (20) is 


Xmd — -D ep /DJ)[ D . 


( 22 ) 


Once x is established, bp can then be determined from 
an MD computation of shear viscosity using Eq. (21) 


(^p)md — 0.53x md 


'1 + 4.76(xmd^-U -1 


(23) 

Here and are the self-diffusion and shear vis¬ 

cosity coefficients from MD simulations. These were com¬ 
puted using the code and techniques described i n PH] and 
[T51 . More information is also provided in Sec.fv] 

Figure |5l shows the results of Eq. (22) (black circles) 
and Eq. ( |23[ ) (red squares). One apparent feature is that 
Xmd begins to differ from 1 as T approaches 1. It is 
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nearly constant over the range 1 < T < 30, then rapidly 
increases for T > 30. This is approximately the coupling 
strength at which there is a known crossover to liquid-like 


behavior DUE]. This is discussed further in Sec. [VI] 


The MD results can be used to check that physical 
requirements are not violated in the modified Enskog 
approach. If the particle size is being modeled as an 
“effective” hard sphere, one should check that the maxi¬ 
mum hard-sphere packing fraction is not exceeded. This 
is known to be fj max = 7t/(3\/2) = 0.74, which gives 
(bp) max = 47r/(3%/2) = 3.0 UU. Indeed, this is ap¬ 
proached, but not exceeded, for the range of data avail¬ 
able in Fig. [5] It is interesting to note that the maxi¬ 
mum physical packing fraction is being approached close 
to the coupling strength at which Wigner crystallization 
is known to occur (r m = 175). The limit of the data 
range is set here by the maximum T value at which we 
can obtain a numerical solution to the HNC equations. 

The MD simulations can also be used to check the 
virial expansion relating bp and \. This was obtained by 
substituting (bp) md from Eq. ( p3| ) for the co-volume in 
Eq. ( [To] ) . Figure [5] shows that this is a reasonable approx¬ 
imation (within tens of percent accuracy) for T < 30, but 
is substantially inaccurate for larger coupling strengths. 
This may be expected from the fact that the virial ex¬ 
pansion is a low-density expansion, which breaks down 
at sufficiently large coupling strength. 


B. A Heuristic Approach 

As discussed before, in Enskog’s theory the probability 
of finding a particle center within a sphere of radius a of 
another particle is zero. Here, we seek an analogy be¬ 
tween the hard sphere system (Fig|2[i) and the Coulomb 



Coupling Strength, T 


FIG. 6. (color online) Co-Volume (bp) and corresponding col¬ 
lision probably enhancement factor (y) computed, (a) using 
the method of Sec. |IV B for three values of a and (b) using 
the methods of Sec. 



Coupling Strength, T 


FIG. 5. (color online) Solutions of Eqs. (221 and (231 obtained 
using MD simulations of self-diffusivity and shear viscosity of 
the OCP. Also shown are the maximum physical co-volume of 
hard spheres (dashed line) and the virial expansion obtained 
using the (&p)md data in Eq. (|10[) (solid line). 


system (Fig[2[)) that provides an effective co-volume in 
terms of coupling strength, bp(T), as well as the collision 
frequency enhancement factor y(r). Although the prob¬ 
ability of finding other particles at any given separation 
does not completely vanish in this situation, the strong 
Coulomb repulsion at close distances creates a region es¬ 
sentially devoid of other particles [18] . 

To quantify this analogy, consider the pair distribu¬ 
tion function g(r ). Physically, ng(r) represents the radial 
density profile around an individual particle. Figure [2]t 
shows example profiles for hard spheres for three values of 
the packing fraction computed using the Percus-Yevick 
approximation m ■ This is a common approximation 
known to accurately represent the pair distribution for 
the hard sphere potential. The figure shows that the 
density is zero within a radial distance of a from the 
particle center. The pair distribution function contains 
the information that particle centers must be at least a 
distance a apart because hard spheres cannot overlap. 

Figure [2[> shows example profiles of the OCP pair dis- 
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FIG. 7. (color online) Co-volume and collision enhance¬ 
ment factors computed from the MD me thod from Sec. |IV A| 
(squares), the heuristic method from Sec. IVB using a = 0.87 


(circles) and from the Stroud-Ashcroft method from Sec. IV C 
(lines). 


tribution function computed from the hypernetted chain 
(HNC) approximation, as well as MD simulations. Like 
the Percus-Yevick approximation for hard spheres, HNC 
is a well established approximation known to accurately 
represent g(r) for Coulomb and screened Coulomb sys¬ 
tems m At weak coupling the density profile ng(r) 
increases over a broad distance, indicating that a defi¬ 
nite exclusion zone is difficult to identify. However, as 
the coupling strength increases into the strongly coupled 
regime, the density profile steepens substantially. Here, 
the exclusion radius can be identified as the radial loca¬ 
tion where g(r) steps from a small value to a value of 
unity order. 

The heuristic approach that we propose is to identify 
the location where g(r) reaches some critical value (a) as 
the effective particle diameter. The co-volume can then 
be computed from this. However, the OCP model has a 
uniform neutralizing background, so the co-volume asso¬ 
ciated with this diameter actually contains two particles. 
Thus, the individual particle co-volume is approximated 
as 

bp~^Trna 3 = ^{a/a) 3 (24) 

where a is defined from g(r = a) = a for some chosen 
value of a. Here, we assume that a is constant (indepen¬ 
dent of r). 

For an actual hard sphere, any value in the range 
0 < a < 1 will provide the same diameter (cr); see Fig. UK 
For a plasma, different values of a will generally give dif¬ 
ferent predictions for the diameter. Here, we determine 
the value of a as that which best represents the OCP 
data under the constraint that it be independent of F. 
Although this value will be chosen by comparison with 
MD data for the OCP, we hypothesize that it represents 


an intrinsic property of Coulomb holes for any given g(r). 
Thus, the value of a is not a fitting parameter. Rather, 
it is presumed to be a universal value that applies to 
ion-ion interactions in any plasma. Although the g(r ) 
profiles will change in different systems, the critical den¬ 
sity described by a will remain the same. This hypothe¬ 
sis will be corroborated by comparison with the Yukawa 
OCP model at different k values in Sec. IV Al Future work 
will test this hypothesis for more realistic systems. Af¬ 
ter bp is obtained in this manner, the virial expansion 
from Eq. (10) is used to estimate \■ This last step is 
not well justified theoretically since the expansion is for 
the hard sphere system, but it is this step that enables a 
“mapping” between the hard sphere model and a plasma. 


Figure [6] shows the co-volume (bp) and collision prob¬ 
ability enhancement factor (x) that result from apply¬ 
ing this method using the HNC approximation for g(r ) 
and three different values for a. Figure [2] shows that at 
weak coupling the radial density distribution gradually 
increases from 0 to 1, so an effective exclusion radius 
is not a well-defined concept in this regime. However, 
Fig. [6] shows that at weak coupling the co-volume is very 
small anyway, so although the expected error in exclu¬ 
sion radius might be large it has a negligible affect on 
the transport quantities in this regime. As the coupling 
strength increases into the moderately coupled regime, 
T > 0.1, the radial density step becomes steeper and 
the co-volume becomes finite. This leads to values of % 
that give rise to order unity corrections to the transport 
coefficients. 


Figure [7] shows that the heuristic method agrees very 
well with the MD evaluation of \ for a ll F < 30 when the 
value a = 0.87 is chosen. The agreement for co-volume 
is less accurate, but is of the order expected from the 
virial expansion (see Fig. [5]). The heuristic method is 
entirely unable to capture the regime T > 30. This may 
demonstrate the limitations of a dense gas picture, since 
it is known that liquid-like behaviors such as caging set 


in this regime m- This is discussed further in Sec. [V| 


C. Relation to Other Works 


Previous work has also considered a correspondence 
with hard sphere systems to gain insight into plasma be¬ 
havior. Ross and Seale m sought a mapping between 
the OCP Coulomb coupling parameter (r) and the hard 
sphere packing fraction (fj) using a variational method. 
The basis of this approach is the Gibbs-Bogolyubov in¬ 
equality, which relates the Helmholtz free energy of a 
given system to that of a reference system, which in this 
case is the hard sphere gas. This was used to approxi¬ 
mate thermodynamics quantities with some success, es¬ 
pecially at large coupling strengths (P > 100) for the 
OCP. A few different results based on this method have 
been published. Stroud and Ashcroft J2D] obtained the 














9 



FIG. 8. (color online) Normalized self-diffusion coefficient of 
the OCP computed from MD simulations (circles), the effec¬ 
tive potential theory based on the Boltzmann-like collision 
operator (triangles), and the modified Enskog method from 
Sec. IV B using a = 0.87 (squares). 


relation 


r = 2fj 1/3 


2 - tj (1 4- 2rj) 2 
2 + fj (1 — fj) 5 


(25) 


This analysis used the Percus-Yevick g(r) to obtain the 
internal energy, and the Carnahan-Starling expression [9] 
to obtain the hard sphere entropy. DeWitt and Rosen- 
feld [21 j provided an alternative expression 


r = 


(1 + 2 fj) 2 
( 1 - 


(26) 


which uses the Percus-Yevick g(r) to obtain both the 
internal energy and hard sphere entropy. More recent 
results by Faussurier and Murillo [22] provide the ex¬ 
pression 


ln(ry) = 31nT — ln(8)+ 


(27) 


[1.845 + 0.006 In P] In 


exp(—1.503 In T +0.22) 

1 + exp(—1.503 In T +0.22) 


which uses an improved Percus-Yevick approximation for 
the internal energy along with the Carnahan-Starling ex¬ 
pression for the hard sphere entropy. 


The results of the co-volume obtained from Eqs. (251 
(27) are shown in Fig. 6b. The x factor obtained from 
applying these in Eq. (101 are also shown. The results 


of this method lead to a much larger estimated packing 
fraction at strong coupling than the method of Sec. |IVB| 
Stroud and Ashcroft have shown that the variational up¬ 
per bound of the excess free energy associated with this 
approach agrees well with Monte Carlo simulation results 
only at very strong coupling (r > 100); see Fig. 2 of [20] . 
Thus, it is likely that the heuristic method underesti¬ 
mates the packing fraction at very strong coupling, but 



Coupling Strength, T 


FIG. 9. (color online) Ratio of the OCP self diffusion co¬ 
efficients calculated using MD simulations and the effective 
potential theory based on the Boltzmann-like collision opera¬ 
tor (diamonds) and including the Enskog correction based on 
the method of Sec. |IV B| for a = 0.5, 0.87 and 0.95. 


the variational method overestimates it for T < 100. Fig¬ 
ure [7] shows that these theories do not accurately capture 
the MD evaluation of the Enskog theory parameters over 
the range of coupling strengths were the effective poten¬ 
tial approach applies. 


V. TRANSPORT COEFFICIENTS 


A. Self-Diffusion 


Figure [8] shows a comparison between MD simula¬ 
tions and the effective potential theory for the OCP 
self-diffusion coefficient. The MD simulation data was 
computed from the Green-Kubo relation as explained 
in m • The theory was computed two ways. The first 
method evaluated the usual effective potential theory as 
described in [I], This applies the potential of mean force 
computed from the HNC approximation to a Boltzmann 
kinetic equation and the associated diffusion coefficient 
from the Chapman-Enskog solution. This corresponds to 
D| P in Eq. (201, and includes up to the second order in 
the Chapman-Enskog expansion; see Eq. (14) of |T] for 
the formula. The second method included the Enskog 
correction, Dp P /\ [see Eq. (20)]. Here, the factor \ was 
computed using the method of Sec. |IV B| with the value 
a = 0.87. 


Including Enskog’s collision probability enhancement 
factor (y) provides a substantial improvement to the ef¬ 
fective potential theory. A subset of the data points 
shown in the figure are given in table [T] along with the 
% difference calculated as — ^ep/xI/^me)- +i§“ 

ure [9] shows the ratio of the theoretical and MD results. 
Without the Enskog correction, the effective potential 
theory overestimates the self-diffusion coefficient by ap¬ 
proximately 30-40% in the range 1 < T < 30. Enskog’s 
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FIG. 10. (color online) Self-diffusion coefficient for the 
Yukawa OOP at four different screening parameters calculated 
using MD (circles), the effecitive potential theory (dashed red 
lines), and the modified Enskog method from Sec. IVB using 
a = 0.87 (solid blue lines). 


FIG. 11. (color online) Normalized shear viscosity of the 
OCP computed from MD: total (circles), kinetic contribution 
(downward triangles) and potential contribution (upward tri¬ 
angles). Three theoretical evaluations are shown: r^ P , t/ep/x 
and Eq. ( |21| |, for which the heuristic method with a = 0.87 
was used to determine x and bp. 


TABLE I. Values of the self-diffusion coefficient computed 
from MD simulations and Eq. (201 where the second order 
Chapman Enskog expansion is used to calculated P . 


r 

Dmd 

d* 2 

% diff 

r 

Dmd 

D* 2 

% diff 

0.075 

246 

237 

3.4 

2.9 

0.482 

0.464 

3.8 

0.1 

133 

131 

1.5 

4.0 

0.337 

0.324 

3.9 

0.2 

31.3 

32.5 

3.7 

5.0 

0.264 

0.255 

3.4 

0.25 

21.6 

21.0 

2.6 

10 

0.131 

0.125 

3.9 

0.3 

14.5 

14.8 

2.1 

20 

0.0653 

0.0622 

4.8 

0.5 

6.02 

5.81 

3.5 

30 

0.0411 

0.0414 

0.6 

0.7 

3.56 

3.28 

7.8 

40 

0.0275 

0.0303 

10 

0.75 

3.16 

2.93 

7.2 

50 

0.0202 

0.0241 

19 

0.9 

2.41 

2.20 

8.7 

55 

0.0180 

0.0215 

20 

1.0 

2.02 

1.87 

7.3 

60 

0.0156 

0.0201 

30 

1.54 

1.07 

1.02 

5.1 

80 

0.00944 0.0151 

60 

2.0 

0.770 0.72 

6.1 

100 

0.00672 

0.0117 

74 


X factor essentially corrects this error if the effective ex¬ 
clusion radius is chosen from an appropriate value of a. 
Here a = 0.87 provides good agreement over the coupling 
parameter range of interest. Other values of a also lead 
to an improved theory, but a value near 0.87 provides the 
most accurate fit. 

The a factor is obtained directly from g(r) based on 
the concept that there is a critical density that defines 
the particle radius. Thus, the numerical value is expected 
to be system-independent. This notion is supported by 
Fig. [lOj which shows the self-diffusion coefficient for the 
Yukawa OCP model at four different screening param¬ 
eters ft = 1,2,3 and 4. Here, the modified Enskog ap¬ 


proach from Sec. IV B is shown to provide a similar ac¬ 
curacy for the Yukawa OCP as the OCP when the same 
value a = 0.87 is chosen. 


B. Shear Viscosity 


Figure [TT] shows a comparison between MD simulations 
and the EP theory for shear viscosity of the OCP. Here, 
three different methods for evaluating the EP theory are 
shown. One is the previous evaluation from |14) based 
on the Boltzmann collision operator ( 77 ^ P ). Here the co¬ 
efficient associated with the first order Chapman-Enskog 
expansion is shown (the formula is given after Eq. (21) 
above). The second is a simple modification that treats 
the lowest-order increased scattering probability factor 
(Vep/x), but not the potential terms in Enskog’s theory 
that arise from the non-local feature of the collision oper¬ 
ator. The third evaluation includes the complete Enskog 
expression from Eq. (21). In addition, the partial con¬ 


tributions from kinetic and potential terms of viscosity 
are shown individually from the MD data; details are 
provided in [14] . 

The figure shows that, once again, application of En¬ 
skog’s increased scattering probability factor (y) im¬ 
proves the effective potential theory over the range 1 < 
r < 30. Applying this correction leads to excellent 
agreement between the theory and kinetic component 
of the MD simulations over the entire range of coupling 
strengths shown. However, the full modified Enskog so¬ 
lution, including the nonlocal effects arising from the bp 
terms in Eq. (21), significantly degrades the accuracy of 


the theory. Furthermore, the method appears to be inca- 
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pable of capturing the minimum in viscosity. The mini¬ 
mum occurs where the potential components exceed the 
kinetic components. These effects are meant to be mod¬ 
eled by the terms in square brackets in Eq. (21), however 
the large T at which this transition occurs coincides with 
where the heuristic approach for obtaining the particle 
size breaks down. 


VI. DISCUSSION OF RESULTS 

The results of sections [IV] and [V] substantiate the sup¬ 
position that a modified Enskog theory can be used to im¬ 
prove the effective potential transport theory, but there 
are limitations. There is particularly strong evidence 
that the increased collision frequency (y) associated with 
the reduced volume that particles of a finite size can 
occupy leads to a significant improvement in the range 
1 < T < 30. Although there is some arbitrariness in how 
the effective particle size is defined (the a factor) even 
a simple heuristic approach essentially removed any dis¬ 
crepancy between the theoretical and MD calculations in 
this regime. It is also supportive that any physically rea¬ 
sonable value of a led to an improvement of the theory, 
although some values were found to lead to better results 
than others. In a future publication, we will test the va¬ 
lidity of the present theory against molecular dynamics 
simulations of realistic plasmas using the techniques re¬ 
cently presented in (2J 

Although the modified Enskog theory led to a substan¬ 
tial improvement in the range 1 T < 30, it essentially 
failed at modeling the potential contributions to viscos¬ 
ity that arise for T > 30. One might suggest that this 
is due to an inaccurate model for the effective particle 
size in this regime. However, if this is the explanation, 
Fig- E shows that the effective particle size must make 
an abrupt adjustment to a much steeper scaling with T 
to extend Enskog’s arguments to this regime. It is hard 
to imagine on a physical basis why the effective particle 
diameter should increase so abruptly at T ~ 30 if the 


medium is to remain dense gas-like. Instead, this may be 
indicative of the breakdown of the dense gas-like picture 
and the onset of a liquid-like regime. 

The data presented here gives further evidence of the 
following analogy between transport regimes of the OCP 
and neutral fluids: (i) Dilute gas like regime (T < 1). 
This is supported by the finding that the dilute-gas Boltz¬ 
mann equation, using an effective interaction potential, 
accurately predicts transport coefficients in this regime, 
(ii) Dense gas like regime (1 < T < 30). Here a Boltz¬ 
mann equation with effective interaction potential gives 
reasonable results, but these are improved by also in¬ 
cluding dense gas effects. In particular, the finite ef¬ 
fective size of particles gives rise to a nonnegligible in¬ 
crease in the collision frequency, (iii) Liquid like regime 
(30 < T < 175). Although there is no abrupt gas-liquid 
phase transition for the OCP, as there is in neutral flu¬ 
ids, there is substantial evidence for liquid-like behavior 
in this regime. Previous work has shown that particle 
caging dominates transport, and that liquid-state scal¬ 
ing relationships, such as Stokes-Einstein, the Arrhenius 
law of viscosity, and other excess-entropy scaling rela¬ 
tionships hold in this regime mug. Additionally, the 
present work has shown that an Enskog-type dense gas 
theory abruptly runs into difficulties as this regime is ap¬ 
proached. (iv) Crystalline regime (r > 175). The OCP 
undergoes a well-known phase transition to a crystal lat¬ 
tice near T = 175. 
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